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We describe a Monte Carlo procedure which allows sampling of the disjoint configuration spaces 
associated with crystalline and fluid phases, within a single simulation. The method utilises biased 
sampling techniques to enhance the probabilities of gateway states (in each phase) which are such 
that a global switch (to the other phase) can be implemented. Equilibrium freezing-point parameters 
can be determined directly; statistical uncertainties prescribed transparently; and finite-size effects 
quantified systematically. The method is potentially quite general; we apply it to the freezing of 
hard spheres. 
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Freezing is the archetypal phase transition, one of the 
prime examples of thermodynamics in action, and a topic 
of ongoing interest ^^ ■ It is therefore remarkable that 
the challenge it presents to computational science has yet 
to be satisfactorily met. The generic problem is to com- 
pute the location of the freezing transition (more gen- 
erally the liquid-solid coexistence curve) on the basis of 
a particle-level model. The approach to this problem 
has evolved little since the pioneering work of Hoover 
and Ree B. The free energy of each phase (fluid (F) 
and crystalline solid (CS)) is computed for states of a 
range of densities, using integration methods which con- 
nect their thermodynamic properties with those of effec- 
tively single-particle reference states, whose free energies 
are known a priori; the two branches of the free-energy 
are then matched to determine the freezing parameters. 
This approach has several drawbacks. The integration 
path may encounter singularities -both real and artificial 
. Corrections may be needed to allow for the fact that 
the path does not quite reach the idealised reference state 
lH . The implicit perspective adopted (that there are two 
separate calculations to be done -one for each phase) has 
meant that predictions for freezing parameters are often 
a synthesis of work done by different authors on different 
system-sizes, making it hard to quantify finite-size effects 

This paper describes a different approach to the prob- 
lem. We build on recent work M], in which we showed 
that the disjoint configuration spaces associated with 
two phases of a many-body system can both be visited 
in a single Monte Carlo (MC) simulation, by harness- 
ing extended-sampling (ES) methods ||l3] to facilitate a 
direct switch from one phase to the other: instead of 
traversing a region where both phases coexist [O] the 
method may be thought of as leaping from one space to 
the other; the role of ES is to allow the system to find the 
'gateway' states from which a leap will be accepted. The 
method was developed |9| to tackle the problem posed by 
two crystalline phases, where interfacial states are com- 
putationally problematic. The same is true of the CS-F 
problem. But significant extensions of the framework 
are needed to address the qualitatively different charac- 
ters of the two configuration spaces. First, the communal 
entropy of the fluid p2| provides a conceptually different 
form of barrier that has to be negotiated to reach the 
gateway states: we show how one can do this. Second, 
the distinct contributing configurations have to be identi- 
fied with care: in so doing we unearth a small but signifi- 
cant flaw, inherent (we think) in all previous simulation- 
studies of CS-phase free energies. The method we de- 
velop is general; we illustrate it here with a study of the 
entropically-driven freezing of hard spheres, where earlier 
studies provide useful benchmarks 0,H,O] 

Consider N particles (hard spheres) confined to volume 
V, variable under a constant external effective pressure 
p [Q, and subject to periodic boundary conditions. The 



configurational weight of a phase may be written as 
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where 7 (CS or F) labels the phase, while 
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Here E is the hard sphere configurational energy WM . 
The prefactor Cq is chosen according to whether the par- 
ticles are regarded as 'strictly classical' (Co = 1) or 'clas- 
sical but indistinguishable' (Cq — -^y The results for 
observables are independent of this choice. The 7-label 
on the integral stands for some configurational constraint 
that picks out configurations {f} that 'belong' to phase 
7. We choose to formulate that constraint as follows 
p5[ . Let Rj . . ■ R'li = {R}^ denote some representative 
configuration of phase 7. Then the constraint may be 
regarded as picking out those configurations which can 
be reached from {R}'^ on the simulation timescale [ [l6[ . 
It is convenient to use the sites defined by {R}'^ as the 
origins of the particle coordinates. Thus we define a set 
of displacement vectors {u} by u^ = r^ — R'J and write 
E''{{u}) = E{{R-' + u}). 

In the case of the F-phase all contributing configura- 
tions are reachable from any one; we may write simply 
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where {R\^ is some specific but arbitrary fluid config- 
uration, which can be selected at random in the course 
of MC exploration of the fluid phase. It is natural to 
choose {R}'~^^ to define the sites of a lattice of the ap- 
propriate symmetry (here fee) and scale flTJl . But one 
must recognise that the complete CS configuration space 
actually comprises a number of distinct mutually inacces- 
sible fragments [Q corresponding essentially to the dif- 
ferent permutations of particles amongst lattice sites [^ . 
By symmetry each fragment should contribute equally to 
the configurational weight; but MC simulation will visit 
(and thus count) only the states within the fragment in 
which it is initiated. The total configurational weight of 
the CS phase is given by multiplying the contribution of 
one fragment by the number of fragments. Since global 
translation (permitted by the boundary conditions) en- 
sures that one fragment includes all possible locations 
of any chosen particle, the number of fragments is the 
number of ways of assigning the 'other' N — 1 particles 
to A^ — 1 Wigner-Seitz cells of some underlying notional 
fixed lattice. This number is not A^! but {N — 1)!. Thus 
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The ratio of the configurational weights of the two phases 
(the ratio of their total a priori probabihties) follows by 
combining Eqs. |l|, ^ and Q 
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from which the Gibbs free-energy-density difference fol- 
lows as 

A.9 = gcs{N,p) - gF{N,p) = - ln7^F,cs (6) 

In writing Eq. M we have chosen to split the fragment 
number {N — 1)! into separate factors of 1/N\ and 1/A^. 
If one so wishes [^ one may regard the former as the 
familiar indistinguishability overcount-correction appro- 
priate for phases (fluids) of non-localised particles. But 
then one must recognise the existence of an analogous 
correction (the 1/N) for the CS phase, in which particles 
are localised — but only relative to one another. It seems 
that this correction has been missed by other authors; 
we shall see that it contributes significantly to finite-size 
effects. 

The relative stability of the two phases is determined 
by the ratio of the associated configurational weights 
(Eq. ^) . To determine that ratio we need a MC procedure 
which visits both solid and fluid regions of configuration 
space. Since, by construction, the system may be trans- 
formed between the CS and F reference states simply by 
switching the representative vectors {R[ ^ Rf^Vi), by 
continuity, any CS (F) configuration 'sufficiently close' 
to the representative one will also transform to a F (CS) 
state under this operation. This phase switch can itself 
be realised as a MC step, so that the phase label 7 be- 
comes a stochastic variable. The set of configurations 
for which the MC switch will be accepted will, however, 
constitute only a small fraction of the respective configu- 
ration spaces. To ensure effective two-phase sampling the 
MC procedure must be biased Q to enhance the prob- 
abilities with which these 'gateway' regions are visited. 
To that end we define an order parameter 



M = M^{{u}) - y {0,[l ~ e{u, ~ u,)] + W{i 



=)} 



Here 9 is the step function. Ti = aui measures the length 
of a notional tether connecting site i to its associated 
particle pO| . Oi measures the overlap (between particle 
i and its neighbours) which would be created by a phase 
switch. The parameter a controls the relative impor- 
tance of Ti and Oi] u^ controls the tether- length domain 
in which each contributes pit] . The equilibrium states of 
both phases are characterised by large M values. The 



'overlap' term contributes in both phases: swapping the 
{R\ vectors will (in general) produce a configuration of 
the 'other' phase in which spheres overlap. The 'tether' 
term contributes only in the F-phase 123] where particles 
may drift arbitrarily far from the sites with which they 
are nominally associated; the tethers provide the means 
to 'pull' the fluid up the communal entropy barrier. We 
identify the gateway states as those which have M = (ie 
Oi = and Ui < Uc, Vi). The constraint that M = im- 
poses on the overlap simply recognises that MC-switches 
which generate overlaps will necessarily be rejected. The 
constraint {ui < Uc) on the tether length is needed to en- 
sure that switches from the fluid create only crystalline 
solid (not defective, glassy) configurations. The entire 
region of configuration space relevant to the problem can 
then be sampled in the multicanonical ensemble defined 
by 

ZiN,p,{v}^T.l dvY[Jdu,e~^-'(i^^'^^ (7) 



D! 



where 

Wiiu}, V) = Emu}) +pV + rj^{M) - S^.cs In (N 

while {rj} represents weights (defined on the M- 
macrostates) which have to be constructed so as to en- 
hance, appropriately, the probabilities of the M = gate- 
way states |24|. Simulation in this ensemble allows one 
to measure the multicanonical probability distribution 
P{M,V,j\N,p,{ri}) from which (unfolding the bias due 
to the weights) one may infer the true equilibrium dis- 
tribution P{M, V, 'j\N,p). The desired ratio of the phase 
probabilities (Eq. ||) follows by 'marginalising' AI and V 
to give the a priori probabilities of the phases. Having 
the underlying distribution of M and V allows one to de- 
termine, in addition, the value of TZp cs ^^ neighbouring 
pressures, using histogram reweighting techniques ] p5[ ]. 

We turn to the MC procedure required for efficient ex- 
ploration of the space spanned by the configuration vari- 
ables {u}, y and 7. It comprises four types of configura- 
tion update, each of which is accepted with a probability 
defined by a Metropolis rule ]g] and reflects the associated 
change in the effective energy H. The first two -particle 
position updates |26|] and volume updates (implemented 
as dilations)- are effected in standard ways Q . The third 
-like the first two- also preserves the phase label; but it 
is novel. In this process, we choose two sites at random (i 
and j say) and identify the corresponding displacement 
vectors Ui and ilj . The candidate configuration is defined 
by the replacements 
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This process can be thought of as an association up- 
date: the particle initially associated with ('tethered to') 
site j is subsequently associated with site i (and vice 



versa) . It changes the representation of the configuration 
(the coordinates {u}); but it leaves the physical config- 
uration invariant. It is required in the fluid phase only 
p7[ . In the fluid phase the particles diffuse far from the 
sites with which they are initially associated; the mem- 
bers of {u} become large and the tethers correspondingly 
so; association updates allow the tethers to respond effi- 
ciently to the influence of the tether contribution to {77}. 
Finally, the 'phase update' (the switch) entails replac- 
ing one set of representative vectors, {R}'^ say, by the 
other, {R}'^ , with the volumes scaled appropriately and 
the displacement coordinates {u} held fixed Pq |. 

Simulations have been performed using systems of TV = 
32, 108 and 256 particles. Figure |] shows the density 
distribution for the N — 256 system in the vicinity of 
the coexistence pressure. Coexistence (Ag — 0; Eq. ||) is 
identified by the equality of the contributions associated 
with each phase (essentially the area under each peak). 

Figure || shows the coexistence pressure for our three 
system sizes plotted as a function of 1/iV [^. The re- 
sults for N = 108 and N = 256 were obtained in the 
fashion just described; the associated uncertainties a[p] 
follow simply from Eq. ^ as a[p\ = a[TZ]/{N \ Av |) where 
Av = [Vp — Vcs]/N and (t[}Z] is the uncertainty in the 
measured ratio of the peak-weights, which is controlled, 
at heart, by the statistics of the inter-phase switch. The 
result for N — 32 was determined differently: this system 
is sufficiently small that transitions back and forth be- 
tween F and CS phases occur spontaneously, over a range 
of pressures, and a density distribution (sampling both 
phases) can be determined -and a coexistence pressure 
inferred -without multicanonical weighting. The three 
points are consistent with the presumed scaling form p9[ . 
The extrapolated prediction (p = 11.49(9)) is, within er- 
ror, in accord with |4| and [|l3[ (see Fig. g inset). 

The lower set of data points shown in figure || gives 
the values of the coexistence pressure implied by our mea- 
surements for A'' = 108 and N = 256 if one fails to fold in 
the 1/N correction in Eq. H. The associated overestimate 
of the CS-configurational weight lowers the predicted co- 
existence pressure by an amount ([In A'^]/[A'^Aw]) which 
is significant for systems of this size, and leads to val- 
ues which it is hard to reconcile (cf the dashed line in 
Fig . H) with the independent measurement at iV = 32 
|P0| . While this correction vanishes in the N —^ 00 limit, 
its existence is potentially important for any systematic 
study of the finite-size scaling of free energies |^ . 

We summarise. We have presented a method which 
allows one to locate liquid-solid coexistence parameters 
(and uncertainties) directly and transparently (Fig. |l|) 
within a single simulation, conducted in the appropriate 
(constant pressure) ensemble. The method avoids the 
need to appeal to integration through to 'distant' refer- 
ence states, double-tangent-constructions or off-the-shelf 
equations of state. It prescribes finite-size effects explic- 



itly and handles systems sufficiently large (cf ||]) that the 
limiting thermodynamic behaviour can be identified with 
some confidence. The method can be readily generalised 
to systems with real (soft) potentials g] and arbitrary 
geometries H. It can also be naturally combined with 
histogram reweighting techniques p3] to allow the full 
coexistence-curve to be mapped efficiently. 
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FIG. 1. The distribution of the density of the system of 
N = 256 particles at pressures (a) just below, (b) at and (c) 
just above coexistence for this N. The mean single phase 
density averages are pp = 0.934(3) and pes = 1.031(4) in 
accord with the coexistence parameters reported in fLSl. 
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FIG. 2. The coexistence pressure for systems of different N 
using Eq. 0both with (•) and without (o) the 1/A prefactor in 
the CS configurational weight. The solid lines is a fit to the 
former; the dashed line is lower by In A/ [A At;]. The inset 
compares our extrapolated value with the results of others, 
with error bars shifted for clarity. 



